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ABSTRACT 

We present, for the first time, a clear A/-body realization of the strong mass segregation solution for the stellar 
distribution around a massive black hole. We compare our A/-body results with those obtained by solving the 
orbit-averaged Fokker-Planck (FP) equation in energy space. The A/-body segregation is slightly stronger than 
in the FP solution, but both confirm the robustness of the regime of strong segregation when the number fraction 
of heavy stars is a (realistically) small fraction of the total population. In view of recent observations revealing 
a dearth of giant stars in the sub-parsec region of the Milky Way, we show that the time scales associated with 
cusp re-growth are not longer than (0.1-0.25) x TriJfh). These time scales are shorter than a Hubble time for 
black holes masses M, <4x 1O 6 M and we conclude that quasi-steady, mass segregated, stellar cusps may be 
common around MBHs in this mass range. Since EMRI rates scale as M~ a , with a £ [1, 1], a good fraction of 
these events should originate from strongly segregated stellar cusps. 

Subject headings: black hole physics — galaxies: nuclei — stellar dynamics — gravitational waves 



1, INTRODUCTION 

The distribution of stars around a massive black hole 
(he nceforth MBH) is a classical problem in stellar d ynam- 
ics dBahcall & Wolfl[l976t [Lightman & Shapirolll977l) . The 
observational demonstration of the existence of nuclear stel- 
lar clusters (henceforth NSCs) — as revealed by a clear up- 
turn in central surface brigthness — in the centers of galax- 
ies makes its study ever more timely. A number of NSCs 
in coe xistence with a central M BH have recently been de- 
tected dGraham & Spitleril 2009) suggesting that NSCs around 
MBHs, like the one in the center of the Milky Way, may be 
quite common. 

The renewed interest in this theoretical problem is thus mo- 
tivated by the observational data in NSCs and, in particular, 
the very rich and detailed data available for the stars orbiting 
the Galactic MBH. At the same time, the prospects for detec- 
tion of gravitational waves (GWs) from extreme mass ratio 
inspirals (henceforth EMRIs) with future GW detectors such 
as the Laser Interferometer Space Antenna (LISA) also urge 
us to build a solid theoretical understanding of sub-parsec 
structure of galactic nuclei. In fact, EMRI rates will depend 
strongly on the stellar density of compact remnants as well as 
on the detailed physics within (9(0.01pc) of the hole, which is 
the region fr om which these inspira lling sources are expected 
to originate dHopman & Alexanderl2005l) . 

Bahcall & Wolfl (ll976l) have shown, through a kinetic treat- 
ment that, in the case all stars are of the same mass, this quasi- 
steady distribution takes the form of power laws, p(r) ~ r" 7 , 
in physical space and f(E) ~ E p in energy space (7 = 7/4 
and p = 7 — 3/2 = 1 /4). This is the so-called zero-flow solu- 
tion f or which the net flu x of s tars in energy spa c e is pre cisely 
zero. iPreto etalJ (12004b and iBaumgardt et all (l2004al) were 
the first to report A/-body realizations of this solution, thereby 
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validating the assumptions inherent to the Fokker-Planck (FP) 
approximation — namely, that scattering is dominated by un- 
correlated, 2-body encounters and, in particular, dense stellar 
cusps populated with stars of the same mass are robust against 
ejection of stars from the cusp. 

The properties of stellar systems that display a range of stel- 
lar masses are only very poorly reproduced by single mass 
models. It is well known from stellar dynamical theory that 
when several masses are present there is mass segregation — 
a process by which the heavy stars accumula te near th e 
center while the lighter ones float outward (Spitzer 1987). 
Accordingly, stars with different mass get distributed with 
different density profiles. By a ssuming a stella r popu la- 
tion with two mass components, Bahcall & Wolfl (11977b — 
hencefort BW77 — generalized their early cusp solution and 
argued heuristically for a scaling relation pi = mi/niH x pn 
that depends on the star's mass ratio only. However, they ob- 
tained no general result on the inner slope of the heavy ob- 
jects; nor did they discuss the dependence of the result on 
the component's number fr action s. On the other hand, it was 
shown long ago by iHenonl (119691) that the presence of a mass 
spectrum leads to an increased rate of stellar ejections from 
the core of a globular cluster, but he did not include the pres- 
ence of a MBH at the center. Henon's work raises the ques- 
tion as to whether multi-mass stellar cusps, obtained from the 
solution of the FP equation, are robust against ejection of 
stars from the cusp. Ejections — due to strong encounters — 
are a priori excluded from the FP evolution, even though they 
could occur in a real nucleus. Furthermore, even if cusps were 
shown by A/-body results to be robust against stellar ejections 
(and we show that they are in this Letter), BW77 scaling can- 
not be valid for arbitrary number fract ions. 

Recently Alex ander & Hopmanl ( 120091) — henceforth 
AH09 — stressed this latter point and have shown via FP 
calculations that, indeed, in the limit where the number 
fraction of heavy stars is realistically small, a new solution 
that they coined strong mass segregation obtains with density 
scaling as Ph(j) ~ r~ a , where a > 2. They have shown 
that there are two branches for the solution parametrized by 
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^ = WW ' T,+M H jM L ■ 'weak branch, for A > 1 corresponds 
to the scaling relations found by BW77; while the strong 
branch, for A < 1, generalizes the BW77 solution. There is 
a straightforward physical interpretation. In the limit where 
heavy stars are very scarce, they barely interact with each 
other and instead sink to the center due to dynamical friction 
against the sea of light stars. Therefore, a quasi-steady state 
forms in which the heavy star's current is not nearly zero 
and thus the BW77 solution does not hold. As A increases, 
self-scattering of heavies becomes important and the resulting 
quasi-steady state forms with a nearly zero current for stars 
of all masses, so BW77 solution is recovered. 

For all these reasons, it is fundamental to verify the 
Bahcall-Wolf solution — as well as its Alexander-Hopman 
generalization — with jV-body integrations. There has been 
a surprisingly small number of j V-body studies o f multi- 
mass systems aro und a MBH (Baumga rdt etal.1 [2004b; 
Freitag et al. 2006), and none of them reported the occurence 
of strong mass segregation. 

In this Letter we use direct A^-body integrations to show 
for the first time that: (i) strong mass segregation is a robust 
outcome of the growth of stellar cusps around a MBH when 
A < 1; (ii) BW77 solution is recovered when A > 1; (iii) 
as a corollary, we conclude that the rate of stellar ejections 
from the cusp is too low to destroy the high density cusps 
around MBHs — even though ejections from the cusp do oc- 
cur. Furthermore, having validated the FP formalism, we pro- 
ceed to use it to estimate the time scales for cusp re-gr owth 
starting from a wider range of models. iMerrittl (120091) ob- 
tained, for Milky Way type nucleus, times in large excess of a 
Hubble time. We use a FP formalism which, in contrast with 
that of the latter author, is tailored to follow the simultane- 
ous evolution of the cusp of different stellar masses without 
any restrictions with respect to the values of f(E) or p(r). 
With our FP solutions we show that, for M. < 5 x 10 6 M Q , 
the times for re-growing stellar cusps are shorter than a Hub- 
ble time. Our results clearly suggest that strongly segregated 
stellar cusps around MBHs in this mass range may be quite 
common in NSCs and should be taken into account when es- 
timating EMRI event rates. 



2. MODELS AND INITIAL CONDITIONS 

We have performed the jV-b ody simulation s with a modi- 
fied version of NBODY4 (lAarsethl U999. 2003) adapted to the 
GRAPE - 6 special-purpose hardware. The code was mod- 
ified to add the capture of stars by the MBH: stars that en- 
ter a critical radius r cap from the hole are captured and their 
mass is added to the hole. The new position and velocity 
of the new massive particle are calculated by imposing that 
the capture process conserves total linear momentum. The 
maximum number of particles supported by the memory of a 
micro-GRAPE board is ~ 1.2 X 10 5 , which have been shown 
to be sufficient to accurately describe the evolution of the bulk 
properties (densities in physical and phase space) of the nu- 
clear stellar cluster dPreto et al.l2004l) . but is not enough to re- 
solve its loss cone dynamics accurately. Therefore, we do not 
attempt a detailed modelling of tidal disruption processes and 
set the capture radius to be equal for all particles. The MBH 
dominates the dynamics inside its influence radius defined 
to be the radius which encloses twice of its mass at t = 0. The 
stellar distribution evolves and reaches its asymptotic quasi- 
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TABLE 1 

N-BODY INTEGRATIONS . I s ' COLUMN: NUMBER OF RUNS ; 2 nd COLUMN: 
SLOPE OF THE DEHNEN'S MODEL INNER CUSP AT t = 0; 3 rd : RATIO OF 
BH MASS TO TOTAL CLUSTER MASS IN STARS; 4 th : f = Nn/N FRACTION 
OF HEAVY MASS PARTICLES; 5 th : ALEXANDER & HOPMAN PARAMETER; 
6 th : INFLUENCE RADIUS r h ; 7 lh : COULOMB LOGARITHM AT r h . THE 
TOTAL NUMBER OF PARTICLES IS N = 1.24 X 10 5 IN ALL RUNS; THE 
MASS RATIO BETWEEN HEAVY AND LIGHT COMPONENTS IS R = 10 FOR 
ALL RUNS. THE TIDAL CAPTURE RADIUS r, 
USE UNITS G = M mc 



mp _ 10~ 7 IN ALL RUNS. WE 
1, WHERE M„„r IS THE TOTAL MASS OF THE 



NUCLEAR CLUSTER AND a IS THE DEHNEN MODEL'S SCALE LENGTH. 

steady state over relaxation time scales (Spitzer 1987): 



T rlx (r h ) ~ 0.34 



CP-phm* In A ' 



(1) 



where 07, and pi, are, respectively, the ID velocity dis- 
persion and sp atial density evaluated at r%. Following 
iPreto et all 42004), we define the Coulomb logarithm In A = 
ln(r h al/2Gm*). 

A realistic mass population of stars with a continuous range 
of stellar masses can be approximately represented by two 
(well-separated) mass scales: one in the range O(IMq) cor- 
responding to low mass main-sequence stars, white dwarfs 
(WDs) and neutron stars (NSs); another with O(10M Q ) rep- 
resenting stellar black holes (SBHs). The relative abundance 
of objects in these mass ranges is overwhelmingly dominated 
by the lighter stars — typical number Tactions of SBHs be- 
ing of O(10" 3 ) dAlexanderi 120051) . The initial stella r model 
is bui lt from a Dehnen model of a spherical galaxy dDehnenl 
[19931) to which a massive particle is added at the center at 
rest dTremaine et ai] 1 19941) . The positions and velocities are 
Monte-Carlo realizations that accurately reproduce the spa- 
tial p(r) and phase space f(E) densities with stars of the same 
mass. In order to generate a two-component model, we pro- 
ceed as follows: (1) we specify the mass ratio R = inn/mi 
between heavy and light stars, and respective number frac- 
tions Jh = Nh/N and f L = 1 -/#, through which the AH09 
A parameter is fixed; (2) we assign the mass m# or nti ran- 
domly to each star according to the statistical weights /# 
and fi, respectively. The resulting model is almost in dy- 
namical equilibrium; deviations of virial ratio from unity are 
< 1 -2%. On a dynamical time scale, phase mixing occurs 
and the virial ratio converges to unity to within a fraction of 
a percent. Following this prescription, the two-component 
models start without any mass segregation, as would be ex- 
pected from a violently relaxed system. Dehnen model's den- 
sity has p(r) oc r" 7 at the center and the corresponding distri- 
bution function f(E) is isotropic. Table 1 gives the list of runs 
and adopted parameters. 

3. FOKKER-PLANCK MODELS FOR SEVERAL STELLAR MASSES 

We also study the evolution of the NSC with a multi- 
mass Fokker-Planck formalism and compare results with the 
A^-body integrations. The time-dependent, orbit-averaged, 
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FIG . 1 . — Evolution of stellar mass within 0.1 r;, from the MBH, for light (3 
upper panels) and heavy (3 lower panels) components. Noisy curves are from 
A'-body integrations; smooth curves are from the Fokker-Planck evolution. 

isotropic, Fokker-Planc k equation in energy space is defined, 
for e ach component ( Spitzei] 119871 : [Chernoff & Weinbergl 
1990), by 



rj7 ,dfi 9Fei 
P ^~dt = ~~dE~ ,E ' i 



'Deej-^-De/i, 



(2) 



D EEJ = 4ir 2 G 2 mltflnA 
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J - q(E) / dE'fj{E') 



IM 

dE'q(E')fj(E') 
f'j 



D E.i = -Y,[^) J E dE'p(E%(E'). 



(3) 
(4) 



In this equation, i,j run from 1 to N c (the number of mass 
components), [ij = m.j/m* where m* = l/N is a reference 

mass. p(E) = 4 f*" m( ® dr r 2 ^2(E-<5>(r)) = -dq/dE is the 
phase-space accessible to each (bound) star of specific energy 
E = -v 2 /2+$(r) > dSpitzerl 19871) . and the total gravitational 
potential <£>(r) is the sum of the contribution from the nuclear 
cluster plus the hole. During our simulations, the stellar dis- 
tribution and its resulting gravitational potential change sub- 
stantially inside r/, only — the region over which the MBH's 
potential is dominant — so we keep the contribution from the 
stars to total $(r) fixed throughout. This system of FP equa- 
tions treats self-consistently both dynamical friction and two- 
body scattering between all components, without any further 



FIG. 2. — The mass density profiles around the MBHs for both components 
at the end of the integrations. A'-body and Fokker-Planck curves are super- 
imposed for comparison. Left panels show ptir) for the light component; 
right panels show pn(r) for heavy component. The arrows signal the location 
of O.lr/, and r/, radii. These plots highlight the asymptotic solution of both 
methods is in good agreement and pn ~ r^ 1H , where yg decreases from val- 
ues > 2 down to ~ 7/4 while moving from the strong to the weak branch of 
the solution. 

approximations other than those inherent to the FP formalism. 
In contrast with Merritt (2009), our treatment is not limited 
to early evolution where the heavy component is just a small 
perturbation on the (time evolving and dominant in number) 
light component. As a result, we can follow both weak and 
strong branches of the solution throughout wi thout restriction. 
We so lve the FP equations (f2l [3) using the Chan g & Cooper! 
( 119701) integration scheme. 

4. RESULTS 

The density of stars around the MBH increases as the cusp 
grows for both light and heavy components until it reaches a 
quasi-steady state; afterwards, the lights start to slowly ex- 
pand. This can be seen in Figure [l] where the mass in- 
side a sphere of O.lr/,, centered on the MBH, is depicted 
as a function of time. This distance is measured with re- 
spect to the instantaneous position of the MBH. Both curves 
from FP and NB are shown together for three different runs 
with A = 0.08,0.23 corresponding to strong segregation and 
A = 13.2 to weak segregation branches. The time scaling 
between NB and FP is related through T^f = InA/N T%f, 
and no further adjustements were made. In the three cases 
shown (as in all others cases tested but not shown), the agree- 
ment between both methods is very good, although there is 
a noticeable tendency for the heavy particles in the NB runs 
to segregate more strongly in the central cusp — this is espe- 
cially the case in the strong branch. Figure Q] also suggests 
that a quasi-steady state (and maximum central concentra- 
tion) have been reached by the end of the runs correspond- 
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FIG. 3. — Evolution of the phase space density (upper panels) and spa- 
tial density (lower panels) for both components. The models starts with 
7=1 /2, purported to represent a cored nuclei; the fraction of stellar black 
holes is / = 10~ 3 and their mass is 10 times larger than that of the light 
stars. Both densities increase monotonically with time and are shown at 
t/T Hx =0,0.05,0.1,0.2,0.25. 



ing to t ~ (0.1 — 0.2)T r i x (rh). We stress that mass segregation, 
whether in the weak or strong branch, speeds up cusp growth 
by factors r anging from 4 to 10 in comparison with the single- 
mass case dPreto et alj2004l) . 

Figure[2]displays the spatial density profiles pi(r) and /?//(r) 
at late times, t ~ 0.27)./ A (r/,). The agreement between both 
methods is again quite good although there is the tendency, in 
the strong branch, for NB's asymptotic slope 7^ to be slightly 
smaller than in FP — for which jL,min = 1.5. The slopes of the 
inner density profiles of the heavy component decrease as the 
solution evolves from the strong to the weak branch when A 
is increased, as expected. In the limit of A >> 1, 7// tends to 
evolve to a quasi-steady state close to the 7/4 solution, while 
for A << 1, 7# > 2. The asymptotic inner density slopes, 
in both solution branches, of the light component extend out 
to ~ O.lr/,; in contrast, the heavy component shows a differ- 
ent behavior depending on the solution branch: on the weak 
branch, 7#'s asymptotic slope also extends only up to ~ O.lr;,, 
while on the strong branch it extends virtually all the way to 
77,. In the strong branch, the density of the heavy component 
exceeds that of the light for r < O.Olr;, (and will therefore 
dominate the interaction events with the MBH); in the strong 
branch, pu > pi throughout. 

Although there are some differences in quantitative detail, 
these NB results broadly confirm the FP calculations and val- 
idate its inherent assumptions — at least in what concerns the 
description of the bulk properties of stellar distributions. A 
more detailed study and description of the stellar dynamics 
around a MBH under different initial conditions and different 
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FIG. 4. — The points represent the relaxation time at the influence radius 
ri, for single mass, cored models of NSCs as a function of MBH mass. The 
shaded area covers the region [0.ir r /. r ,0.27" r / v ], which is time for cusp re- 
growth computed with the Fokker-Planck equation. Its finite width results 
from the distribution of parameters: the number fraction of stellar BHs / £ 
[10~ 3 , 10" 2 ] and mass ratio R = 10, 15, 20. The horizontal curve (with points) 
corresponds to a Hubble time 13 Gyr. NSCs with MBH masses below that of 
SgrA* can re-grow their stellar cusps after < 6 Gyr. 

models (allowing for more than the two components reported 
here), stellar ejections (which do occur in the NB runs) and 
captures, and further comparisons between NB and FP meth- 
ods, are outside the scope of this Letter and are the subject of 
another work in preparation. 

5. IMPLICATION FOR GALACTIC NUCLEI AND SOURCES OF GWS 

The analysis of the number counts of spectroscopically 
identi fied, old stars in the sub-parsec region of our own Milky 
Way dBuchholz et al.1 120091; iDo etal.ll2009l) — believed to be 
complete down to magnitude K = 15.5 — reveals a deficit of 
old stars with respect to the high number a strongly segregated 
cusp would entail. Although the slope of the density profile is 
still weakly constrained, the best fits from num ber counts data 
seem to exclude with certainty slopes 7 > 1 (Schod eTet al.l 
2009), and there could be a core with a stellar density de- 
creasing towards the center, 7 < 0, although such a fit is only 
marginally better than one with 7 ~ 1/2. 

Although we deem to be too early to conclude for the in- 
existence of a segregated cusp around SgrA*, since the de- 
tectable stars (essentially giants) are still a small fraction of 
the stellar population as a whole, we next compute the time 
necessary for cusp growth if at some point a central core is 
carved in the stellar distribution. Having validated the FP ap- 
proach and its results we study equations (0 |3), which are 
orders of magnitude faster to solve than NB integrations. 

We choose as initial condition a model with 7=1/2 (which 
is the minimum slope that allows an isotropic solution around 
a MBH), since the isotropization time — the time necessary for 
the establishment of this shallow cusp s tarting from a hole in 
the spatial distribution — is <C T r ix(Xh) (Merritt 2009), and we 
are interested in the evolution over 0{T r i x ) time scales. Fig- 
ure [3] shows the evolving phase-space f(E) and spatial p{r) 
densities for both components (R = 10 and / = 0.001 consti- 
tute our fiducial case). It can be seen that, by t ~ 0.25 T^r;,), 
cusps with 7£ ~ 1.5 and 7// ~ 2 (or pi ~ 0.05 and pn ~ 0.5 
in phase space) are fully developed; a little earlier, at t ~ 
0.2 T r i x (yh), the density cusp Ph(/) is already fully developed 
down to r ~ O.Olr;, (~ 0.02 pc if scaled to the Milky Way 
nucleus). If there was some event carving a hole in the stel- 
lar distribution around SgrA* more than 6 Gyr ago, then there 
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was enough time for a very steep cusp of stellar BHs to have 
re-grown. 

The number fraction / of SBHs is sensitive to the initial 
mass function (IMF) of high mass stars. The re are indica- 
tions the IMF in galactic nuclei is top-heavy dManess et al.l 
2007) so we adopt a range of values / G [10" 3 ,10 -2 ]; the 
mass distribution of SBHs i s also weakly constrained so we 
follow lO'Learv et al.l (120091) in considering several mass ra- 
tios R = 10, 15 and 20. Figure |4] shows the relaxation times 
at the influence radius for nuclei with MBH masses in the 
range of interest for LISA; the straight line is a linear fit 
to the points. The shaded region corresponds to the range 
[0.1 T r i x {rh)i 0.2 T r i x (ri,)] and represents the time stellar cusps 
take to grow starting from an isotropic core. The shaded re- 
gion's width results from the distribution of values for R and 
/. In the ranges we adopted, increasing R or / both have 
the effect of decreasing the time for cusp growth. At early 
times, SBHs essentially evolve under dynamical friction with 
characteristic time scale Tjf ~ T r i x /R; increasing / leads to an 
increased rate of self-scattering between SBHs at late times. 

6. SUMMARY AND DISCUSSION 

Our results show that strong mass segregation is a robust 
outcome from the growth of stellar cusps around MBHs. 
We have used A^-body integrations with two masses — light 
and heavy components representing main sequence stars and 
stellar BHs respectively — , and compared the results with 
those obtained with the FP formalism. The broad agree- 
ment between both methods validates the FP description of 
the bulk properties of time-evolving stellar distribution around 
a MBH — and its underlying assumptions. The differences of 
quantitative detail are the subject of another work in prepara- 
tion. 

Using the FP equation to study cusp growth under a vari- 
ety of initial conditions purported to represent cored nuclei, 
we have shown that the time scales associated with cusp re- 
growth are clearly shorter than a Hubble time for nuclei with 

3 GRACE: see http://www.ari.uni-heidelberg.de/grace 



MBHs in the mass range M, <5x 1O 6 M — even though the 
relaxation time, as estimated for a single mass stellar dis- 
tribution, exceeds a Hubble time in the upper part of this 
mass range. Therefore, our work strongly suggests that quasi- 
steady — strongly segregated — stellar cusps may be common 
around MBHs with masses in this range. 

EMRIs of compact remnants will be detectable 
by LISA precisely for MBHs in this mass range 
dde Freitas Pacheco et al.l 120061: lAmaro-Seoane et ail 120071 : 
Baba k etaU 120071) . Estimates for event and detection rates 
by LISA cost umarily assume that t he stel lar cusps are in 
steady state dHopman & Alexander! I2006alfb1) . But recent 
observations reveal a dearth of giants inside 1 pc from SgrA* 
and raise the possibility that cored nuclei are common — this 
scenario has been thoroughly explored by lMerritil (120091) . 

Our results strongly suggest that stellar cusps can re-grow 
in less than a Hubble time. The existence of cored nuclei still 
remains plausible though — especially for nuclei with MBHs 
in the upper part of the mass range — , since time scales are 
still quite long (e.g. 6 Gyr in Milky Way type nuclei). How- 
ever, since EMRI rates scale as M~ a , a € [4,1], and re- 
growth times are < 1 Gyr for M. < 1.2 X 10 6 M Q , we still 
expect that a substantial fraction of EMRI events will origi- 
nate from segregated stellar cusps. Finally, indirect observa- 
tions alone will reveal whether there is a "hid den" cusp of old 
stars and their dark remn ants around SgrA* dWeinberg et al.l 
120051 IPreto &~Sahll2009l) . 
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